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Abstract 

A free energy functional for a crystal that contains both the symmetry conserved and symmetry 



U . 

< 

broken parts of the direct pair correlation function has been used to investigate the crystallization 



of fluids in three-dimensions. The symmetry broken part of the direct pair correlation function 
has been calculated using a series in ascending powers of the order parameters and which contains 
three- and higher-bodies direct correlation functions of the isotropic phase. It is shown that a very 
accurate description of freezing transitions for a wide class of potentials is found by considering 
the first two terms of this series. The results found for freezing parameters including structure of 
CJ | the frozen phase for fluids interacting via the inverse power potential u(r) = e (a/r) n for n ranging 

from 4 to oo are in very good agreement with simulation results. It is found that for n > 6.5 the 

> - 

f^ . fluid freezes into a face centred cubic (fee) structure while for n < 6 the body centred cubic (bec) 

m , 

structure is preferred. The fluid-bcc-fcc triple point is found to be at 1/n = 0.158 which is in good 
,— ^ ' agreement with simulation result. 
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I. INTRODUCTION 

Freezing of a fluid into a crystalline solid is a particular, but an important example of a 
first-order phase transition in which the continuous symmetry of the fluid is broken into one 
of the Bravais lattices. The transition in three-dimensions is marked by large discontinuities 
in entropy, density and order parameters; the order parameters being proportional to the 
lattice components of one particle density distribution p(r) (see Eq.(2.3)). Efforts have been 



made for over six decades 



QB 



2J to find a first principle theory which can answer questions 
such as, at what density, pressure and temperature does a particular liquid freeze ? What 
is the change in entropy and the change in density upon freezing ? Which of the Bravais 
lattices emerges at the freezing point for a given system and what are values of the order 
parameters? 

A crystal is a system of extreme inhomogeneities where value of p(r ) shows several orders 
of magnitude difference between its values on the lattice sites and in the interstitial regions. 
The density functional formalism of classical statistical mechanics has been employed to 
develop theories for freezing transitions l2|, |3j. This kind of approach was initiated in 1979 
by Ramakrishnan and Yussouff (RY) [4| which was latter reformulated by Haymet and 
Oxtoby J5j. The central quantity in this formalism is the reduced Helmholtz free energy of 



both the crystal, A[p], and the liquid , A(j>i) [2j. For crystals, A[p] is a unique functional of 
p(r) whereas for liquids, A(pi) is simply a function of liquid density p\ which is a constant, 
independent of position. 

The density functional formalism is used to write an expression for A[p] (or for the 
grand thermodynamic potential) in terms of p(r) and the direct pair correlation function 
(DPCF). Minimization of this expression with respect to p(r) leads to an expression that 
relates p(r) to the DPCF. The DPCF that appears in these equations corresponds to the 
crystal and is functional of p(r) and therefore depends on values of the order parameters. 
In the RY theory the functional dependence of the DPCF on p(r) was neglected and was 
replaced by that of the coexisting liquid of density p\. Attempts to improve the RY theory 
by incorporating a term involving three-body direct correlation function of the coexisting 

n n 

liquid in the expression of A[p] have failed 6|, \A . The efforts made by Tarazona [8|, Curtin, 

121 ] in the direction of developing a theory using 



Ashcraft and Denton 19, 10] and others 



laU 



ii 



what is referred to as the weighted density approximation have also met with limited success 



only. 

The reason, as has been pointed out recently 13j, ll4j], is that at the fluid - solid transition 
the isotropy and the homogeneity of space is spontaneously broken and a qualitatively new 
contribution to the correlation in distribution of particles emerges. This fact has been used 
to write the DPCF of the frozen phase as a sum of two terms; one that preserves the 
continuous symmetry of the liquid and the other that breaks it and vanishes in the liquid. 
An exact expression for the free energy functional was found by performing double functional 
integration in density space of a relation that relates the second functional derivative of A[p] 
with respect to p(r) to the DPCF (see Eq. fl2.7D ). This expression of free energy functional 
contains both the symmetry conserved and the symmetry broken parts of the DPCF. 

The values of the DPCF as well as of the total pair correlation function (described in 
Sec II) in a classical system can be found from solution of integral equation, the Ornstein 
- Zernike (OZ) equation, and a closure relation that relates correlation functions to pair 



potential 15]. The integral equation theory has been quite successful in getting values of 



pair correlation functions of uniform liquids 15j, but its application to find pair correlation 
functions of symmetry broken phases has so far been limited. Recently Mishra and Singh [16j 
have used the OZ equation and the Percus - Yevick (PY) closure relation to obtain both the 
symmetry conserved and symmetry broken parts of pair correlation functions in a nematic 
phase. In the nematic phase the orientational symmetry is broken but the translational 
symmetry of the liquid phase remains intact whereas in a crystal both the orientational 
and the translational symmetries of the li quid phase are broken. Since, closure relations are 



derived assuming translational invariance 15), they are valid in normal liquids as well as in 



nematics but may not in crystals. In view of this, Singh and Singh |l3j suggested a method 
in which the symmetry broken part of the DPCF is expanded in ascending powers of order 
parameters. This series contains three- and higher - bodies direct correlation functions of 
the isotropic phase. The first term of this series was evaluated and used in investigating 
the freezing transitions in two- and three-dimensions of fluids interacting via inverse power 
potentials |13l . |l7J and freezing of hard spheres into crystalline and glassy phases [14| . It has 
been found that contribution made by the symmetry broken part to the grand thermody- 
namic potential at the freezing point increases with softness of the potential |l3L Il7| . This 
suggests that for long - ranged potentials the higher order terms of the series may not be 
negligible and need to be considered. 



In this paper we calculate first and second terms of the series (see Eq.(2.29)) which involve 
three and four-bodies direct correlation functions of the isotropic phase. We calculate the 
four-body direct correlation function by extending the method developed to calculate the 
three-body direct correlation function. The values found for the DPCF are used in the free- 
energy functional and the crystallization of fluids is investigated. We show that all questions 
posed at the beginning of this section are correctly answered for a wide class of potentials. 

The paper is organised as follows: In Sec. II we describe correlation functions in liquids 
and in crystals and calculate them. The symmetry broken part of the DPCF is evaluated 
using first two terms of a series in ascending powers of order parameters. These results 
are used in the free-energy functional in Sec. Ill to calculate the contributions made by 
different parts of the DPCF to the grand thermodynamic potential at the freezing point. 
In Sec. IV we calculate these terms and locate the freezing points for fluids interacting 
via the inverse power potentials and compare our results with those found from computer 
simulations and from approximate free energy functionals. The paper ends with a brief 
summary and perspectives given in Sec. V. 

II. CORRELATION FUNCTIONS 

The equilibrium one particle distribution p(r) defined as 

P<?) = (YJ<?-?i)), (2-1) 



J2 S (r -ri)) 



where r\ is position vector of the I th particle and the angular bracket, (....), represents the 
ensemble average, is a constant, independent of position for a normal liquid but contains 
most of the structural informations of a crystal. For a crystalline solid there exists a discrete 
set of vectors Ri such that, 

P(r) = p(r + %), forall ^ (2.2) 

This set of vectors which appears at the freezing point due to spontaneous breaking of 
continuous symmetry of a liquid, necessarily forms a Bravais lattice. The p(r) in a crystal 



can be written as a sum of two terms: 

p(^)=p + P ( V) (2.3a) 

where 

p®P) = Y, *>&&*■ ( 2 - 3b ) 

G 

Here po is the average density of the crystal and pa are the order parameters (amplitude 
of density waves of wavelength 2n/\G\). The sum in Eq. (12.3bj) is over a complete set of 



reciprocal lattice vectors (RLV) G with the property that e * = 1 for all G and for all 
Ri. We refer the first term of Eq. (l2.3ap as symmetry conserved and the second as symmetry 
broken parts of single particle distribution p(r) . 

The two-particle density distribution p^ 2 '(r 1,^2) which gives probability of finding simul- 
taneously a particle in volume element dri at r\ and a second particle in volume element 
<ir 2 at r2, is defined as 



J2Y,5(?i-?j)H?2-?k))- (2-4) 



p {2) (r 1 ,r 2 ) 
The pair correlation function g(r\, 7^2) is related to p^(r 1,^2) by the relation, 

»<?.,*,) = ^M- (2-5) 

p{r 1 )p(r 2 ) 
The DPCF c(^r 1, 7^2)) which appears in the expression of free-energy functional A[p] is 

related to the total pair correlation function /i(ri, 7*2) = g{r 1^2) — ~\- through the Ornstien 

- Zernike (OZ) equation 



c(r 1 , r 2 ) = h(r 1 ,r 2 ) - / dr z c(r u r 3 )p(t 3 )h(t 2 , r 3 ). (2.6) 

The second functional derivative of A[p] is expressed in terms of c(r"i, 7*2) as [2J 

^H -*(-"-» -c( ?1 ,? 2 ), (27) 



5p(ri) <5p(r 2 ) p(ri) 

where 5 is Dirac function. The first term on the right hand side of this equation corresponds 

to ideal part A-^fp] of the free energy whereas the second term corresponds to excess part 

A ea .[p] arising due to interparticle interactions. 



In a normal liquid all pair correlation functions defined above are simple function of 
number density p and depend only on magnitude of interparticle separation | r 2 — r\\ = r. 
This simplification is due to homogeneity which implies continuous translational symmetry 
and isotropy which implies continuous rotational symmetry. In a crystal which is both 
inhomogeneous and anisotropic, pair correlation functions can be written as a sum of two 
terms; one that preserves the continuous symmetry of the liquid and the other that breaks 
it |ialld|. Thus 

h(?i,? 2 ) = h {0 \\?2 ~ ?i\,Po) + &<Vi, r 2 ; [p]) (2.8) 

c(?i,? 2 ) = c (0) (|^2 - ?i\,Po) + c<Vi,?25 [p])- (2.9) 

While the symmetry conserving part (h^ and c^) depends on the magnitude of inter- 
particle separation r and is a function of average density po, the symmetry broken parts h^ b ' 
and c^ b > are functional of p(r) (indicated by square bracket) and are invariant only under a 
discrete set of translations corresponding to lattice vectors Ri, 

h®{? u t 2 ) = h b (^ 1 + %,t 2 + R^ (2.10) 

c (6) (n, r 2 ) = c 6 (^ + R h t 2 + %) (2.11) 

If one chooses a centre of mass variable r c = (ri + r 2 )/2 and a difference variable 

r^ = ~f 2 — "ri, then one can see from Eqs. (12.101) and ( 12. lip that h> b > and c- b > are periodic 

unctions of the centre of mass variable and a continuous function of the difference variable 



18|. Thus 



h®(? 1 ,? 2 ) = Y t e* d * c ti G >(?) (2.12) 

G 

c W(^ 2 ) = ]Te fe c( G )(^) (2.13) 

G 

Since h^ and c^ are real and symmetric with respect to interchange of r\ and -r 2 ; 
h(~ G ^(r) = h,( G \~r) and h <yG \— 1 r) = hS G '(^) and similar relations holds for c^(r). 

Substitution of values of h(ri,'r 2 ) and c(ri, r 2 ) given by Eqs. ( 12. 8 p and ( 12. 9 p in Eq. 
( 12. 6 p allows us to split the OZ equation into two equations; one that contains h^\ c^ and 
po while the other contains hS b \ c^ and p(r-3) along with h^°\ c^ and po : 



hW(\?2-? 1 \) = cW(\?2-?i\)+f> /^3C (0) (|^3-^l|)/l (0) (l^3-^2|) (2-14) 



and 



fr (6) 0?i, r 2 ) = c<Vi, ? 2 ) + J ^ 3 c (0) (|r 3 - nD^ft) - Po)^ (0) (l^3 - r 2 \) 
+ jd? 3 p(t 3 ) [c^(f x ,f 3 )h^(\f 3 -? a \) + cW(|r 3 - ^D^Hn^s) 
+c(Vi^3)/i (6 Vi l r f 3)]. (2.15) 

Eq. (j2.14p is the well known OZ equation of normal liquids. We use it along with a 
closure relation to calculate the values of these correlation functions and their derivatives 
with respect to density po- The derivatives of c^(r) are used to find values of three- and 
four- bodies direct correlation functions of the isotropic phase. 

Eq. (I2.15P is the OZ equation for symmetry broken part of correlation functions. In 
order to make use of it to find values of h^ and c~ b ' for a given p(r) we need one more 
relation (closure relation) that connects h^ b ' with c- b >. Alternatively, if we know values of 



one of these functions then Eq. (12.151) can be used to find values of the other function 19]. 
Here we calculate c^(r 1, 7*2) using a series in ascending powers of (p(r) — po). 

A. Calculation of h^°\ c-°> and their derivatives with respect to p 

We use the OZ equation (12.141) and a closure relation proposed by Roger and Young 
which mixes the Percus-Yevick (PY) relation and the hypernetted chain (HNC) relation in 
such a way that at r = it reduces to the PY relation and for r — y 00 it reduces to the 
HNC relation and is written as 



h^°'(r) = exp[—/3u(r)] 



exp[x(r)f(r)} 



(2.16) 



f(r) 

where x( r ) = h^°'(r) — c^°'(r) and f(r) = 1 — exp(-ipr) is a mixing function with 
adjustable parameter < if) < 00, to calculate pair correlation functions and their derivatives 
with respect to density p. The value of if) is chosen to guarantee thermodynamic consistency 



between the virial and compressibility routes to the equation of state 20]. 



The differentiation of Eqs. (j2.14p and ( I2.16p with respect to p yields the following rela- 
tions, 



and 



dh^{r) dc®(r) 



dp 



+ /^V°)(r> (0) (l"-r|) 



dp 

+p /^^rl/,(o) ( |^_^) 

J op 

J op 



dh(0) ( r ) r o , m r < \*t m^W 

— exp[— pu{r)\exp[x{r)j{r)\- 



dp 



dp 



(2.17) 



(2.18) 



d 2 h^\r) d 2 c^{r) 



dp 2 



dp 2 



+ 2 / dr ' 



&<0, M, ( ».(l ? '- ? l) + c( o) M ^ (0, (P'-?l) 



and 



+ p I dr' 



dp " dp 

dc^(r')dh^(\f-t\) (0) , d 2 h^(\r'-t\) 
~^~p dp + ° (r) dp 

+ q 2 ^ Hk ~r\) 



(2.19) 



g 2 /lW(7 



exp[— {3u(r)]exp\x(r) f (r)] 



d 2 x{r) fdx(r) 



dp 2 



+ 



dp 



f(r) 



(2.20) 



The solution of the closed set of coupled equations (12.141) and (J2.17p - (j2.20p gives values 



of/t(°)(r), c(°)(r), 
u(r). 



9fe(°)(r) dcWjr) 8 2 fe(")(r) 



rJp 



9p 



V 



and 



— rr-r^ as a function of r for a given potential 



The pair potential taken here are the inverse power potentials, u(r) = e {<j/r) n where e, a 
and n are potential parameters and r is the molecular separation. The parameter n measures 
softness of the potential; n = oo corresponds to hard-sphere and n — 1 to the one component 
plasma. The reason for our choosing these potentials is that the range of potential can be 
varied by changing the value of n and the fact that the equation of state and mel ting curves 
of these potentials have been extensively investigated by computer simulations 2ll428| for 
several values of n so that "exact" results are available for comparison. The more repulsive 



8 



(n > 7) systems have been found to freeze into a face-centred cubic (fee) structure while the 



soft repulsions n < 7 freeze into a body-centred cu 
fee triple point is found to occur at - ~ 0.15 25|, 



:>ic crystal (bec) structure. The fluid-bec- 



26 



. The atomic arrangements in the 



two cubic structures are very different; the fee is close-packed in real space and the density 
inhomogeneity is much sharper than for the bec which is open structure in real space but 
close-packed in Fourier space. However, in spite of this difference in the atomic arrangements, 
the two structures have small difference in free energy (or chemical potential) at the fluid 



solid transition 



25 



- 28j and therefore a correct description of the relative stability of the 
two cubic structures is a stringent test for any theory. 

The inverse power potentials are known to have a simple scaling property according to 
which the reduced thermodynamic properties depend on a single variable which is defined 
as 

7 = p a 3 (/3e) 3/n = p*T* ( ~ 3/n) 

where /3 = l/fc^T; ks is the Boltzmann constant and T temperature. Using the scaling 
relation the potential is written as 

4tt V /3 1 



0<r) = y—ij - 

( \V3 

where r is measured in unit of ao = ( -^- ) 

In Fig. 1 we plot values of c^(r), c » (r ' and — | a for n = 6 and 7 = 2.30 which is 
close to the freezing point. 

B. Calculation of three- and four-body direct correlation functions 

The higher-body direct correlation function is related to the derivatives d m c ( -°\r)/dp m 
as follows I2I: 



^^)=/^ 3 4 0) (n,^,?3;Po), (2.21) 



<9 2 c (0) (r,p ) [ r-> 4 0) (ri,r 2 ,r 3 ;po) 
ar 3 



dp' 



dp 
dr 3 I dr 4 cf\r lj ^ 2 ^3,r 4 ;p ), 



(2.22) 



etc., where c m are m-body direct correlation function of the isotropic phase of density 
Pq. These equations can be solved to find values of Cm by writing them as a product of pair 
functions. For c 3 (^1,^21^3) one can write as J6[, 



.(o) 



c 3 ; (n, r 2 , r 3 ) = t(r 12 )t(r 13 )t(r 2 3), 




(2.23) 

where a line linking particles i and j denotes a t(r) function and each circle ( representing 
a particle) carry weight unity. The value of t(r) is found from the relation f)2.2ip . 



dc(°>(r,p ) 
dp 




(2.24) 



where the half -black circle represents the particle over which integration is performed over 
its all configurations and all circles carry weight unity. Using known values of dc^ (r, p )/dpo 
we solve this equation to find values of t(r) for different density po{or 7) following a method 
outlined in ref 6[|. The values of t(r) as a function of r are shown in Fig. [2] for n = 6, 4 and 
7 = 2.30, 5.60, respectively . 

Taking derivative of both sides of Eq. (I2.23P with respect to p one gets, 



&SVi,?„?.) __ st^) ^ + i(ri2) ^r H ) i(rB) + i(ri2)((ri3) s*M. (2 . 25) 



<9p <9p 

Substitution of this in Eq. (12.22]) leads to 



<9p 



dpo 



d 2 c (0) ( 



dpo 2 



dr 



dt{r) t{r')t{\r - ?|) + t(r)^^t(|^ - ?|) + £(r)£(r') ^'^ " '' ' ' 



dpo 



dpo 



dpo 



(2.26) 



10 



where r 12 = r, r 13 = r' and r 23 = \~r — ~r\. As values of t(r) are known, Eg. (12.261 ) is used 
to find values of dt(r)/dpo in same way as Eg. ( 12.241) was used to find values of t(r). In Fig. 
[3] we plot dt(r)/dp for n = 4, 6 and 7 = 5.60, 2.30. 

Guided by the relation of Eq. ()2.24p we write dt(r)/dp as 



dp 



s(r) / <ir s(r")s(|r — r|), 



-o 

a 



(2.27) 



where a dashed line connecting particles z and j is s(r) function. Using the already 
determined values of dt(r)/dpo at a given value of po (or 7) we determine values of s(r) in 
same way as values of t(r) were determined from known values of dc^(r)/dpo. In Fig. H]we 
plot values of s(r) for n = 6, 4 and 7 = 2.30, 5.60 as a function of r. 

From Eqs.([2]22D, (12T25J) and ( lOTj) we get 



3 4 

9 ©_....-_ .0 Q - 

,(°)f 



C4 (ri,r 2 ,r 3 , r 4 ) 



b if b 6 



21 21 



(2.28) 



where a dashed line represents s(r )-bond and a full line t(r)-bond. We calculate values 
of c 3 and c A and plot them in Appendix. 



C. Evaluation of c^(r\, r^) 

The function c^(r^i, 7^2) can be expanded in ascending powers of (p(r) — po) as J2J, Il3| . 

c (fe) (ri,r 2 ; [p]) = / dr 3 cf ) (j 1 ,r 2 ,r 3 ;p )(p(j 3 ) - p Q ) 

1 f _> /" _> (0). 



+ 2 l T3 dV 4 ^ ri ' r2 ' r3 ' r4 ' /? °) (jC '( r3 ) ~ /°o)(p( 7 "4) - Po) 

+ ... , (2.29a) 



11 



(2.29b) 



where black circles represent integration over all configurations of these particles and each 
carries weight p(ri) — Po = *YucPG e%G ' r% whereas each white circle carries weight unity. In 
writing Eq.f l2.29bj) use has been made of Eqs.flZSD and (ESED- 

Usefulness of series of Eq.(2.29) depends on how fast it converges and on our ability 
of finding values of c m . We have already described the calculation of c 3 and c\ . The 
same procedure can be used to find c m for m > 4. We, however, find that for a wide 
range of potentials it is enough to consider the first two terms of the series (2.29). In 
fact, for most potentials representing the inter-particle interactions in real systems one may 
need to consider the first term only as contribution made by the second term to the grand 
thermodynamic potential at the freezing point turns out to be negligibly small unless the 
potential has a long range tail. 

1. Evaluation of first term of Eq. (2.29) 

Substituting value of {p{r^) — po) from Eq.(2.3) and using notations, r = ri — r\, 
r = r% — n, r c = 5(7*1 + r 2J we find 

= c (M) (ri, r 2 ) = ^p G e i5 -" c t(r)e-5^^ f dr't(r')t(\r' - r |)e^'. (2.30) 

This is solved to give [13j, ll4| 



c (M) (ri,^ 2 ) = £e iG ^£cpV)Wr)^(G), (2.31) 




G Im 



where 



(G,l)/ \ 

c, r 



0oZ)I>i(Ma,O.7J, (l Gr ) B h(r,G). (2.32) 



h h 



12 



Here ji(x) is the spherical Bessel function, Y[ m (x) the spherical harmonics, 



A 1 (/i,I a ,0 = (») I1+ *'(-l) 



(2f! + l)(2Z 2 + l) 
(21 + 1) 



[C g (h,l 2 ,l;0,0,0)f 



and 



(2.33) 



B h (r,G) = 8t(r) / dkk 2 t(k) Jh (kr) / dr'r'Hir'^kr'^Gr') 



(2.34) 

where C g is the Clebsch-Gordan coefficient. The crystal symmetry dictates that / and 
l\ + I2 are even and for a cubic crystal, m — 0, ±4. 

The values of c\ ' (r) depend on order parameters p^ = po A*G5 where p^ — e_G an d 
on magnitude of G. In Figs. [5], [6] we plot and compare values of c\ ' (r) for bcc and fee 
crystals at the melting point for potential n = 6, 7 S = 2.32, a fecc = 18 and aj cc = 32 (see 
Tabled]). The values given in these figures are for the first and second sets of RLV's. As 
expected, the values are far from negligible and differ considerably for the two structures. 
The value is found to decrease rapidly as the value of / is increased; the maximum contri- 
bution comes from I = 0. We also find, as shown in Fig. [7] , the value of c\ ' (r) decreases 
rapidly as the magnitude of G vector increases; the maximum contribution comes from the 
first two sets of RLV's. The other point to be noted is that at a given point r, values of 
c\ ' (r) are positive for some G vectors while for others the values are negative leading to 
mutual cancellation in a quantity where summation over G is involved. 



2. Evaluation of second term of Eq. (2.29) 

The contribution arising from the second term of Eq.(2.29) is sum of three diagrams in 
which the last two contributions are equal. Thus, 



: (6 ' 2) (ri,^ 2 ) = i 




(2.35) 



If we write r = r4 — r"i and 7^4 = r" +^r c — \r and use other notations defined above, the 
first diagram can be written as 



13 



4 3 



c^'Vi,^) 



^(OEE^^ 31 



+G 2 )-(rc-±r) 



G\ G2 



dr't(r')t(\r' - r\)e lGl - r / dr" s{r")s{\r" - r\)e 



(2.36) 



This is solved to give 



^(n,^ 



E^ e EEffi(** 



iml^J, 



G Im I'm' 



where 



(2.37) 



Hml'm'V ) 



E ** A* E E A Sk™ 2 M h (r, Gi)M fa (r, tf) 

Gi limi ^2»Tt2 



Jl'l^J^mxCGl)^^^). 



(2.38) 



Here K = G-G U 



Ji'hk 



A"'1'2 = lfi > (?' 

Izmz 



S^ ( n\h+h+l' (^Y 



(2/ 1 + l)(2/ 2 + l)(2/ / + l) 
(2/ + 1) 



1/2 



C 9 (/i, Z 2 , Z 3 ; 0, 0, 0)C g (Z', Z 3 , 1; 0, 0, 0)C fl (Zi, Z 2 , k; m 1 , m 2 , m 3 )C g (l', Z 3 , Z; m', m 3 , m); 



(2.39) 



M h {r,GJ = / dr , r' 2 j, 1 (Gr , )t(r') / dkk 2 t(k)j h (kr)j h (kr') (2.40) 



and 



M, 2 (r,X) = / dr"r" 2 j h (Kr")s(r") / dkk 2 s(k)j h (kr) 3l2 (kr") . (2.41) 



The crystal symmetry dictates that all Zj are even and for a cubic crystal all to* are and 



±4. 



From the second diagram of Eq.(2.35) we get 
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4 3 



& b 



EC^'Vl,^ 



^)EE^p G2 ^ 



+G 2 ).(^ c -|r) 



Gi G2 



dr,s(r f )t(\r-^\)e tG ^ r d?"s(r")s(\?" -?'\)e iG ^ r . (2.42) 



This is solved to give 



c ^'Vi^ 2 ) = E" G 



e « . c ^ ^ c lml , m , (r)Y lljn , (G)Yi 



lm(T), 



G Im I'm 1 



where 



(2.43) 



(G,2,2) , n 



Gl h"ll /2"12 £3™-3 



ii'la^l^mi^)^™.^ ( 2 - 44 ) 



j^ll'hhh __ QOf? 

mm'mim2m3 V 



Jl+^2+i'/ 1\i' 



(2/ 1 + l)(2/ 2 + l)(2/' + l) 



1/2 



(2/ + 1) 
C<?(^i, ^2, ^3; 0, 0, 0)C g (l', l 3 , 1; 0, 0, 0)C g (h, l 2 , h] m 1 , m 2 , m 3 )C g (l', l 3 , I; m', m 3 , m); 

(2.45) 



N hM3 (r,G,G 1 ) = t(r) / dr'r >2 S (r% 1 (G 1 r')B h (r , ,K)A l3 (r,r'), (2.46) 



and 



A h (r,r') = / dkk 2 t(k)ji 3 (kr)j h (kr') 



(2.47) 



B h (r',K) = / dr"r" 2 j h (Kr")s(r") / dkk 2 s(k)j h (kr% 2 (kr"). (2.48) 

The total contribution arising from the second term of Eq.(2.29) is 
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c (6 ' 2) (n,^ 2 ) = E e ^EE [^3(0 + cSW] *UW m ,(G), (2.49) 

G Zm Z'm' 

where /, /' are even and m = 0, ±4 for cubic lattices. In Figs. [8] and |9] we plot values of 

(G,2) / \ _ (G,2,l)/ x (G,2,2)/ x / 9 r n x 

C W'm' l > r J ~~ C lml'm> V ) + C lml'm' V)l l/.OUJ 

as a function of r for bcc and fee structures for n = 6, 7 S = 2.32, afe cc = 18 and a/ cc = 32. 
The values given in these figures are for the first two sets of RLV's for I — V — and 2, and 
m = m' = 0. These are the terms which mostly contribute to c Zm ^,(r); the contributions 
from terms / ^ V and m ^ m' are approximately order of magnitude smaller. For a bcc 
lattice we find two sets of values; one for G vectors lying in the x-y plane and the other for 
the rest of the vectors. Since all vectors of the first set of RLV's of a fee lattice are out of x-y 
plane we get only one set of values. For the second set of RLV's of a fee lattice, though two 
sets of values are found but they are close unlike the case of bcc lattice where the two sets 
of values differ not only in magnitude but also in sign. The values differ considerably for the 
two cubic structures. The value of c im ^, (r) decreases for both bcc and fee structures rapidly 

^ ( C 1 1 'l 

as the magnitude of G vectors increases as was found in the case of c\ ' (r). Furthermore, 
the values of c Zm j,^,(r) at a given value of r is positive for some G vector and negative for 
others. 

In order to compare magnitude of contributions made by the first and second terms of 
Eq.(2.29) we calculate c^ G ' l \k^6 kl (f) k ) and c( G ' 2 \k,0 k ,(f)k) defined as 

cW(k,e k ,<f> k ) = pJ2 J drcf' l \r)e^Y lm (r)YL(G) 

lm 

= 4tt pJ2 ^ l Y lm (k)Y; m (G) / dr r 2 c^ l \r) 3l {kr) (2.51) 

Jo 



lm 



and 



lm I'm' 

= ^ PE E (0'>W*)>7m'(3) / dr r 2 c^l'(r)Ji(kr). (2.52) 

lm I'm' ^° 
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In Figs. [TU] and [TT] we compare using colour codes (shown on the right hand side of each 
figure) the values of these functions arising from the first and second terms of Eq.(2.29) for 
both fee and bec structures for n = 6, 7 S = 2.32, af cc = 32 and «b cc = 18. The values 
given in Fig. [10] are for a fee lattice for a G vector of first and a vector of second sets, i.e. 
Gia = 4.25, 9 Gl = 54.7° and <f> Gl = 45° and G 2 a = 4.91, 9 G2 = 0°, <p G2 = 90°. The values 
of ka are taken equal to 4.25, 4.91 which are magnitude of Gia and G 2 a , respectively. 
In Fig. [TT] we compare the values of c^ G '^(k, 9k, 4>k) and c^ G ' 2 \k, 9k, 4>k) for a bec lattice for 
dao = 4.37, 9 Gl = 90°, <f> Gl = 45° and G 2 a = 6.19, 9 G2 = 90°, <p G2 = 0° and ka = 4.37 
and 6.19. From these figures it is clear that the contribution made by the second term to 
c^(ri, r^) is small compared to the first term indicating fast convergence of the series. As is 
shown below, in the expression of free energy functional, c^^i,^) is averaged over density 
and order parameters and also there is summation over G vectors, As a consequence, the 
contribution of second term of Eq.(2.29) is found to be order of magnitude smaller than the 
first term. We show that the consideration of the first two terms of Eq.(2.29) is enough to 
give accurate description of freezing transitions for a wide class of potentials. 

III. FREE-ENERGY FUNCTIONAL AND LIQUID-SOLID TRANSITION 

The reduced free-energy functional A[p] of a symmetry broken phase can be written as 
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iiU 



A[p] = A ld [p]+A^[p]+A^[p] (3.1) 

where 

AM = I d?p(?) [ln{p(?)A) - 1] (3.2) 



A^[p] = A ex {pi) + P(ii ~ HpiA)) / d?(p(?) - Pl ) 



and 



\ I ctfi / rfr 2 (p(^) - pi){p{? 2 ) - pi)c^(\r 2 7, | ) (3 3) 



A®[p] = -\ ld?x fd?2(p(?i) - Po)(p(r 2 ) - po)c (6) (n,^ 2 ) (3.4) 
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Here A is cube of the thermal wavelength associated with a molecule, /3 = (fcgT) -1 , k% 
being the Boltzmann constant and T is the temperature, AiJ(pi) is excess reduced free energy 
of the coexisting isotropic liquid of density pi and chemical potential p and p = p\ (1 + Ap*) 
is the average density of the solid. 

lf X\? 2 -? 1 \) = 2j dXxf d\' c® (\? 2 -^ pt + XX' (po-pi)) (3.5) 

and 



c (6) (n^ 2 ) =4 / dXX f dX' [ dtt f d£c ( » fa^M'poi&PG 
Jo Jo Jo Jo v 



(3.6) 



The expression for the symmetry conserving part of reduced excess free energy A ex [p] 
given by Eq.f l3.3l) is found by performing double functional integration of [13|, Il7| 



5 Aex [p] (0)/|-> -± |\ (o *7\ 

-- -c y> [\f 2 - fi|). (3.7) 



5p(r 1 ) 5p(r 2 ) 

This integration is carried out in the density space taking the coexisting uniform fluid of 
density p\ and chemical potential p as a reference. The expression for the symmetry broken 
part A e x [p] given by Eq. (13.41) is found by performing double functional integration of [13j, [17| 



o Aex[p\ (6)/-+ -*■ \ to o\ 

-- -c (> (r u r 2 ), (3.8) 



5p(rx) Sp(r 2 ) 

in the density space corresponding to the symmetry broken phase. The path of integration 
in this space is characterized by two parameters A and £. These parameters vary from to 
1. The parameter A raises the density from zero to the final value po as it varies from to 
1, whereas parameter £ raises the order parameter from to its final value pa- The result 
is independent of the order of integration. 

In locating the transition the grand thermodynamic potential defined as 

-W = A-(3p J drp{r) (3.9) 

is generally used as it ensures the pressure and chemical potential of both phases remain equal 
at the transition. The transition point is determined by the condition AW = Wi — W = 0, 
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where W, is the grand therrnodynarnic potentia! of the co-existing h qu id. The expression of 

AW is found to be |l3|, 



AW 



dr 



1 



p(r)ln 



Pi 



- (p(r) ~ Pi) 



- dr± dr 2 {p(r 1 ) - pi)(p(r 2 ) - pi)c {0) (\r 2 - »*i|) 

- / dr ! / dr 2 (p(r 1 ) - p )(p(r 2 ) - p )c (6) (n>^2)- 



(3.10) 



Minimization of AW with respect to p(r) subject to the perfect crystal constraint leads 



to 



In 



p{r } 



Pi 

where 



0+ / d? 2 (p(? 2 ) - Pl )(PW2 -Fi|) + / rfr 2 (p(^ 2 ) -p )c w (?i,? 2 ), (3.11) 



i*W/ 



and 



c(°)(|^ 2 -f 1 |)= / dAcW(|^ 2 -^i|;p, + A(po-pi)) 



£ (b) (n,r 2 )= /" dX f dfc® fa, ? 2 ;Apd, fro)) 

Jo Jo 

The value of the Lagrange multiplier <j) in Eq. ( 13.111) is found from the condition 



, tf&± = 1 
VJ po 



(3.12) 



where V is volume of the system. 

It may be noted that, in principle, one needs only values of symmetry conserved and sym- 
metry broken parts of the DPCF to determine p(r) that minimizes the grand potential W. 
In practice, however, it is found convenient to do minimization with respect to an assumed 
form of p{r) . The ideal part is calculated using a form of p(r) which is a superposition of 
normalized Gaussians centred around the lattice sites, 
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p(r) 



a \ 3 / 2 v-^ 

I 2_, ex v 



TT 






—a [~r — R, 



(3.13) 



where a is the variational parameter that characterizes the width of the Gaussian; the 
square root of a is inversely proportional to the width of a peak. It thus measures the non- 
uniformity; a = corresponds to the limit of a uniform liquid and an increasing value of a 
corresponds to increasing localization of particles on their respective lattice sites defined by 
vectors Ri. For the interaction part it is convenient to use the expression of p(r) given by 
Eq.(2.3). The Fourier transform of Eq. fl3.13p leads to pc = PoPg-, where pa — e ~ G //4 ° '• 

A. Evaluation of c(°)(r) and c^' (7*1,7*2) 

The values of c^°\r) for a given liquid density pi and the average crystal density po 
are found from the known values of c^- ' (r, p) where p varies from pi to po by performing 
integrations in Eq.( l3.5p which can be rewritten as 



c (0 \r,p ) 



2 I dXX f d\'c {0) (r; Pl (l + AA' Ap*] 



(3.14) 
'0 Jo v 7 

where Ap* = (p — Pi)/pi- The integrations have been done numerically using a very fine 

grid for variables A and A . Since at the freezing point p/Ap* << 1 one can use Taylor 

expansion to solve Eq. fl3.14j) leading to 



c^(r,p ) = c^(r, Pl ) + \p^P M ^ r pi Pl) + O (p 2 Ap* 2 ) 



(3.15) 



Since the order parameters that appear in c^(r 1,^2) are linear in c^' 1 ' (7*1, 7^2) and 
quadratic in c*- 6 ' 2 ^(r*i, 7*2), the integration over £ variables in Eq. (l3.4j) can be performed 
analytically leading to 



C ( Vl,?2) 



E«' 



G.7 C 



J2^ G ' 1] ^ Y i*m(G)Y lm (r) 



lm 



_(G,2,1) / x _(G,2,2) 



+ EEISM +^C(r))Y; m ,(G)Y lm (r) 

lm I'm' 



(3.16) 
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where 



-(G,l)/ \ 

c, r 



l PG J2J2 A ^ 1 ^ 1 ^ (\ Gr ) B h(r,G), (3.17) 



Q hl2 (r,G,G l )Y* mi (G 1 )Y* m2 (K), (3.18) 



C lm,l'm' V/ r- 

Gi l 1 m 1 Z 2 m 2 



_(G,2,2) / \ 

C lm,l'm'V) n 



£ Z2 PG ^ PK Yl Yl 2Z A mSm 2 m3^' ( lp T 
G\ hmi l2m 2 hm-j, 

N lll2h {r,G,G 1 )Y l l m2 {G 1 )YC zm3 {k) 1 (3.19) 



with 



-l ri 



B h (r, G) = 2 I dXX I d\'B h (r, G; XX' p) , (3.20) 



-1 /•! 



Q llh (r,G 1 G 1 )=2J dXxj d\'Q hh (r,G,G i; AA>) , (3.21) 



Q hh (r, G, Gr, p) = M h (r, G x] p) M h (r, K; p) 



and 



N hl2h (r,G,G 1 ) = 2 J dXxJ dX' N hhh fr,G,G x ;XX' p) • (3.22) 

The quantities Bi^r, G), Mi x (r, G±), Mi 2 (r,K) and A^ 1 / 2 / 3 (r, G, Gi) are defined by Eqs. 
(I2.34p . ( I2.40p . f)2.4ip and (I2.46P respectively. The integrations over A and A have been 
performed numerically by varying them from to 1 on a fine grid and evaluating the functions 
Bi 17 Qi^ 2 and Ni t i 2 i 3 on these densities. Since these functions vary smoothly with density 
and their values have been evaluated at closely spaced values of density the result found for 
cs b >(lri, -r 2 ) is expected to be accurate. 
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B. Evaluation of AW 

Substituting expression of p(r) given by Eqs.(2.3) and (I3.14p and of c^(r) and c^(ri, r*_) 
riven above in Eq. (l3.10p we find 



where 



AW AWm AWo AW { b l} AWl 

N N N N N 



(3.23) 



AW, 



id 



N 



1-(1 + A 7 ) 



+ Inpi - -In ( 



3, fa 

n 



(3.24) 



AW 1 40) 1 2 ^ ,2^<0) 



G^O 



(3.25) 



AW, 



(i) 



TV 



G G 2 ' 



(3.26) 



AWv! 2) 1 



A' 



±p,(l + A 7 ) 2 £ £ W-o-a$* (G2 + -G 



G G 2 



■d,G,2) (-+ i I7* 

2 



(3.27) 



where A7 = (j s — ji) /ji] the subscripts s and Z stand for solid and liquid, respectively. 
Here AW ic [, AW , AW^ and AW^ are respectively, the ideal, symmetry conserving and 

symmetry broken contributions from first and second terms of series (2.29) to AW. The 

_____ ______ — > — > — > 

prime on summation in Eqs. fl3.26p . ( I3.27P indicates the condition G 7^ 0, G\ 7^ 0, G 2 7^ 0, 

G + G 1 ^ and G + G 2 + and 



^ 0) (G)= I d7# \r,>n)c«' ' 



(3.28) 



^ G,1) fe + i3) =^ G ^^A 1 (/ 1 ,/ 2 ,/)y^(G)|^ J , 2 Qo 



(3.29) 
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r» (S, + 13) =*<*■> (3, + i3) + 2* <w > (3 2 + 13 



c(G,2,l) 
C 



G2+ 2 G j - gL^^LLL / , A mm>m,m? Y l*m>( G ) Y l*rm ( G l) Y hm 2 ( K ) 
' G\ Ira I'm' l\m\ l 2 m2 

Jd?j h Qgt) g, 1 , 2 (r,G,G' 1 )e<^ + ^)^ m (f), (3.30) 



C(G,2,2) 

C 



G 2 + ~G J - - 2_^/i Gl /iK 2^ 2^ 2^ 2^ 2^ ^mm4 3 1 m 2 m 3 ^'m'( G ') y i2m2( G 'l)^3m3(- ft ') 
d Zm I'm' lim\ l 2 ra 2 Z3ITI3 

Jdr Jh (±gA Ar y2i3 (r,G,G 1 )e< 32+ ^)-^ m (f). (3.31) 



The terms — ^, — ^— and — jf— are respectively second, third and fourth orders in order 
parameters. 



IV. RESULTS FOR LIQUID-CRYSTAL TRANSITION 

We use above expression of AW/N to locate the liquid - fee crystal and the liquid - bec 
crystal transitions by varying ji , A7 and a. For a given ji and A7, AW/N is minimised 
with respect to a; next A7 is varied untill the lowest value of AW/N at its minimum is 
found. If this lowest value of AW/N at its minimum is not zero, then 7; is varied until 
AW/N = 0. The values of transition parameters, 7/, A7 and a for a given lattice structure 
can also be found from simultaneous solution of equations Q ,^ , (^r) = 0, ^ (^r) — 
and AW/N = 0. 

In Tabled we compare values of different terms of AW/N (see Eq. fl3.23[) ) at the freezing 
point for potentials with n = 4,6, 6.5, 7, 12 and 00. The values corresponding to hard spheres 
are taken from ref. [l4|. The contribution made by the symmetry broken part to the grand 
thermodynamic potential at the freezing point is substantial and its importance increases 
with the softness of the potential. For example, while for n = 00 the contribution of the 
symmetry broken part is about 8% of the contribution made by the symmetry conserved 
part, it increases to 45% for n — 4. As this contribution is negative, it stabilizes the solid 
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phase. Without it the theory strongly overestimates the stability of the fluid phase especially 
for softer potentials. This explains why the Ramakrishanan - Yussouff theory gives good 
results for hard core potentials but fails for potentials that have soft core and/or attractive 
tail. 

The other point to be noted from these results is about the convergence of the series 
(2.29) which has been used to calculate ^'(ri,^)- The contribution made by the second 
term of the series to the grand thermodynamic potential at the freezing point is found to be 
negligible compared to that of the first term for n > 6 and for n < 6, though the contribution 
is small but not negligible. For example, while for n = Q this contribution is about 2% of 
the first term, for n = 4 this increases to 18%. From these results one can conclude that the 
first two terms of the series of Eq.(2.29) are enough to describe the freezing transition for a 
wide class of potentials. 

In Table HIl we compare results of freezing parameters 7/, 7 S , A7, Lindemann parameter 

Pa 3 

L n and where P is the pressure at the transition point, of the present calculation wit 

e 



those found from computer simulations 2ll428l] and with the results found by others 
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30- 



32| using approximate free energy functionals. The Lindemann parameter is defined as the 
ratio of the mean field displacement of a particle to the nearest neighbour distance in the 
crystal. For the fee crystal with the Gaussian density profile of Eq. ( 13.13!) it is given as 



1/2 
2 — ) (4-1) 



a fcc a 



>V3 • 



where a,f cc = (4/p ) is the fee lattice constant. For the bec crystal, 

/ 2 \ 1/2 

L n = -2— > (4-2) 

where a\, cc = (2/p ) is the bec lattice constant. 

In Fig. [12] we plot 72 vs 1/n at the transition found from simulations and from the present 
calculations. 

One may note that simulation results have spread (see Table [XT]) and do not agree within 
each others uncertainties. This may be due to application of different theoretical methods 
used in locating the transition and system sizes in the calculations. The other sources of 
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errors include the existence of an interface, truncation of the potential, free-energy bias, etc. 



Agrawal and Kofke 25| who have reported results for < 1/n < 0.33 have considered a 
system of 500 particles only. Since they have not used finite size corrections, their results 
for softer potentials (say n < 6) may not be accurate. For example, they reported that 
for 1/n > 0.16 fluid freezes into a bcc structure but for 1/n = 0.25 they found that 7/ 
for the fluid - bcc transition is higher than that of the fluid-fcc transition. The recent 
calculations where large systems have been considered 26j-[28| results are available for n > 5 
(or 1/n < 0.2). From these results it is found that fluid freezes into fee crystal for n > 7 
and for n < 7 the bcc structure is preferred; the fluid-bcc-fcc triple point is estimated to be 
close to 1/n ~ 0.15. 

From Table [Til and Fig. [12] we find that our results are in very good agreement with 
simulation results for all cases. We find that for n > 6.5 the fluid freezes into fee structure 
while for n < 6 it freezes into bcc structure. The fluid-bcc-fcc triple point is found at 
- = 0.158 (see the inset in Fig. IT2l) . The value of Lindemann parameter found by us is, 



however, somewhat lower than those found by Agrawal and Kofke [25j and Saija et al [29 ]. 
The energy difference between the two cubic structures at the transition is found to be small 
in agreement with the simulation results [28 1. 

V. SUMMARY AND PERSPECTIVES. 

We used a free energy functional for a crystal proposed by Singh and Singh [13] to 
investigate the crystallization of fluids interacting via power law potentials. This free-energy 
functional was found by performing double functional integration in the density space of a 
relation that relates the second functional derivative of A[p] with respect to p(r) to the 
DPCF of the crystal. The expression found for A[p] is exact and contains both the symmetry 
conserved part of the DPCF, c^\r,p) and the symmetry broken part 0^(7^1,7^2)- The 
symmetry conserved part corresponds to the isotropy and homogeneity of the phase and 
passes smoothly to the frozen phase at the freezing point, whereas the symmetry broken 
part arises due to heterogeneity which sets in at the freezing point and vanishes in the liquid 
phase. The values of c-°\r) and its derivatives with respect to density p as a function of 



interparticle separation r have been determined using an integra 
the OZ equation and the closer relation of Roger and Young 



equation theory comprising 
20j. From the results of c a ^ r ' 
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and — c „ -i we calculated the three- and four-bodies direct correlation functions of the 

dp 

isotropic phase. These results have been used in a series written in ascending powers of the 
order parameters to calculate c^(ri, r^)- The contributions made by the first and second 
terms of the series have been calculated for bcc and fee crystals. The contribution made by 
second term is found to be considerably smaller than the first term indicating that the first 
two terms are enough to give accurate values for c^ b '(jri, r^)- The values of c^(r) for bcc 
and fee structures are found to differ considerably. 

The contribution of symmetry broken part of DPCF to the free energy is found to depend 
on the nature of pair potentials; the contribution increases with softness of potentials. In 
case of power law potentials we found that the contribution to the grand thermodynamic 
potential at the freezing point arising from the second term of the series (2.29) which involves 
four-body direct correlation function is negligible for n > 6 and small but not negligible for 
n < 6. For n = 4 the contribution made by the second term is about 18% of the first 
term. The contribution made by second term is positive whereas the contribution of the 
first term is negative. As the net contribution made by the symmetry broken term is 
negative, it stabilizes the solid phase. Without the inclusion of this term the theory strongly 
overestimates the stability of the fluid 
reported in this paper and elsewhere [141 . 
theory gives good results for hard core potentials but fails for potentials that have soft core 
/or attractive tail. 

The agreement between theory and simulation values of freezing parameters found for 
potentials with n varying from 4 to oo indicates that the free energy functional used here 
with values of c^-r i, r^) calculated from the first two terms of the series (2.29) provides an 
accurate theory for freezing transitions for a wide class of potentials. Since this free energy 
functional takes into account the spontaneous symmetry breaking, it can be used to study 
various phenomena of ordered phases near their melting points. 

Acknowledgements: One of us (A.S.B.) thanks University Grants Commission (New 
Delhi, India) for award of research fellowship. 
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Appendix A 

In this appendix we calculate c 3 (r*i, 7*2,^3) and c\ (7*1, ^"2,^3,^4)- Using the notation 
r = |7*2 — Ti|, r = 1^3 — ?"i| and |r — r| = |r3 — r2| we write c 3 (t"i, r2, T3J as (see 
Eq.(2.23)) 

4 0) (r,^') =t(r)t(r')t(|r / -^|). (Al) 

The function i(|r — r\) can be expanded in spherical harmonics, 



t(\r'-?\) = -J2Mr,r')Y lm (r)Y i ; n (?), (A2) 

Im 

where, 



^•00 
Ai(r,r')= dqq 2 t(q)ii(qr)ji(qr'). (A3) 

Here j;(x) is the spherical Bessel function and Yj m (r) the spherical harmonics. 
From Eqs.(jAT]) and flA2l we get 



4V,/) = -J] A(r,r')>/ m (f)F^(f'), (A4) 

Im 

where 



Di{r,r') = Ai(r,r')t(r)t(r'). 
The Fourier transform of Eq. (1A4l) defined as 



gives 



cfHqxM = 327r^(-i)'A(gi,g 2 )^ m (gi)^(g 2 ), (A5) 
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where, 



A (?i,92) = P dr r I dr' r' ji{qir)ji{q 2 r')Di(r,r') 



(A6) 



,(0) 



The value of c 3 (g 1; q 2 ) is plotted in Fig. [13] for qi = q 2 = q ma x for various angle 9 such 
that < \q x + q 2 \ < 2q max , where 9 is angle between ~q 1 and ~q 2 as shown in the figure. 
The values plotted in this figure correspond to qa>o = 4.3 and for n = 6, 7/ = 2.30 (full line) 
and n = 12, 7^ = 1.17 (dashed line). In Fig. [TJ]we plot values of c 3 Hj} 1 ,~q 2 ) f° r equilateral 
triangle with various side lengths. The values for n = 12, 72 = 1.17 are in good agreement 
with the values given in ref |6| (see Figs. 3 and 4 of ref [6|). 

For c\ (7^1,^2,^3,^4) the contribution arises from three diagrams shown in Eq. (l2.28j ). 
Using the notation | r4 — r 1 1 = r", \r^ — r 2 \ = \r — r\, \r ^ — r^\ = \r — r \ and other 
notations defined above we get 







(A7) 



Each diagram of Eq. (1A7|) has two circles connected by three bonds- two s - bonds (dashed 
line) and one t - bond (full line), one of the remaining circles is connected by two t - bonds 
and the other by two s - bonds. By permuting circles one can convert one diagram into 
another. The values of c\ (r ', 7*, 7* ) depend on three vectors 'r^ 1 and ~r" . 

We calculate c\ C9i, 92, 93) defined as 



4 (9i,9 2 ,9 3 ) = P 



Id? f df' fdf" e-^e-^'e-^'c^^yy), (A8) 



Using Eq. (1A7j) and writing each diagram in terms of t and s bonds we get 



4(91, 9 2 > 9 3 ) 



^ E £ E(- ? ) (/1+ " + " )A ^r 3 ^ 1 ^(9i,9 2 ,9 3 ) 

limi hm2 ^3"i3 
Kms (cO^imi (92)^ 2 m 2 (&) + Vi WO^ms (fe)^^ (&) 

+(-l) h ^ imi (g 1 )y, 3 * m3 (g 2 )y i2m2 (g 3 )] , 



(A9) 
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where 



\m1m2m3 
2 \i 



ihh 



(2h + l)(2l 2 + l) 
4tt(2/ 3 + 1) 



1I/2 



and 



C g (h, h, h; 0, 0, 0)C g (h, hi h\ mi, m 2 , m 3 ), (A10) 



/*oo /»oo /*oo 

Mi 1 i 2 i 3 (qi,q 2 ,q 3 ) = p 3 / dr r 2 s(r) dr' r'H{r') I dr" r" 2 s(r") 

Jo Jo Jo 

3h (li r )Jh (l2r')ji 2 (q 3 r")A h (r, r')E h (r, r"). 
A^ir,^) is defined by Eq. (lA3|) . Ei 2 (r,r") is given as 



(AH) 



Ei{r,r") 



dq q 2 s(q)ji(qr)ji(qr"). 



(A12) 



The values of c\ (^1,^2^3) depend on magnitudes and directions of vectors ~q±, ~q 2 an d 
q 3 . In Figs. [15] and [16] we use color codes(shown at the right hand side of each figure) to 
plot values of c\ (q^ ~q 2 i ~q s ) for q 1 = q 2 = q% = q max as a function of (J) q2 and gs for different 
choices of 6 q2 and Q qz . The values of g max a is taken equal to 4.3 as in Fig. [T3] While the 
values plotted in Fig. [151 correspond to 9 qi = 0°, the values plotted in Fig. [TBI correspond 
to 9 qi = 90° and <p qi = 0°. These figures show how the values of £4 (q 1 ,~q 2 ,~q 3 ) depend 
on orientations of vectors q x , g* 2 an d ~q%- Emergence of ordering in maxima and minima 
depending on orientations of these vectors is evident. 
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TABLE I: Freezing parameters 7; , A7 and the contributions of ideal symmetry conserving and 
symmetry broken parts arising from first and second terms of Eq. (2.29) to AW/N at the transition 
point. 



n 


Lattice 


ll 


A7 


AW id /N 


AWo/N 


AW^/N 


AW { b 2) /N 


4 


bcc 
fee 


5.57 
5.60 


0.007 
0.008 


2.86 
3.52 


-2.09 
-2.64 


-0.94 
-1.03 


0.17 
0.15 


6 


bcc 
fee 


2.30 
2.32 


0.011 
0.012 


2.56 
3.48 


-1.99 
-2.75 


-0.58 
-0.72 


0.01 
0.002 


6.5 


bcc 
fee 


2.04 
2.03 


0.014 
0.013 


2.38 
3.34 


-1.89 
-2.69 


-0.50 
-0.66 


0.001 
0.001 


7 


bcc 
fee 


1.86 
1.84 


0.015 
0.014 


2.29 
3.39 


-1.85 
-2.76 


-0.44 
-0.63 


0.000 
0.000 


12 


fee 


1.17 


0.034 


3.71 


-3.14 


-0.57 


0.000 


00 


fee 


0.937 


0.106 


4.44 


-4.10 


-0.34 


0.000 



TABLE II: Comparison of parameters 7/ , 7 S , A7 , the Lindemann parameter L and the pressure P 
at the coexistence found from different free-energy functional and computer simulations. MWDA 
stands for modified weighted density approximation, RY DFT stands for Ramakrishnan - Yussouff 
Density functional theory. MHNC stands for modified hypernetted-chain closure relation and 
MSMC for Mayer sampling Monte Carlo. 



n 


Lattice 


Theory/ Simulation 


ll 


Is 


A7 


L 


Pa 3 
e 


00 


fee 


Present result 


0.937 


1.036 


0.106 


0.09 


11.46 






MWDA-static reference [TJ] 


0.863 


0.964 


0.115 


0.13 








MWDA [32j 
RYDFT [3ft I3J 


0.906 


1.044 


0.116 


0.10 








0.980 


1.146 


0.174 


0.06 








Simulation [22fl 


0.939 


1.037 


0.104 


~0.13 








Simulation [23J 


0.942 


1.041 


0.105 










MC Simulation [25J 


0.94 


1.041 


0.107 


0.12 


11.70 



Continued on next page 
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XAjjL/E II — Continued from 


previous page 




n 


Lattice 


Theory/ Simulation 


ll 


Is 


A7 


L 


Pa 6 

€ 






MC Simulation 


29 


0.939 


1.037 


0.104 




11.57 


12 


fee 


Present result 


1.17 


1.21 


0.034 


0.11 


23.67 






MWDA-static reference [121] 


1.12 


1.16 


0.037 


0.14 








MWDA/MHNC [3J 
MC Simulation* [25j 


1.19 


1.26 


0.059 


0.10 








1.17 


1.22 


0.042 


0.14 


23.64 






MSMC technique [2J 
MC Simulation |26j 


1.16 


1.20 


0.037 




23.24 






1.16 


1.21 


0.037 




23.41 


7 


fee 


Present result 


1.84 


1.87 


0.014 


0.12 


64.97 






MC Simulation* [2J 


1.85 


1.88 


0.017 


0.15 


64.98 






MC Simulation ]2$ 


1.84 


1.87 


0.016 




64.22 




bee 


Present result 


1.86 


1.89 


0.015 


0.18 


67.12 






MC Simulation [2jj 


1.83 


1.86 


0.015 




63.88 


6.5 


fee 


Present result 


2.03 


2.06 


0.013 


0.12 


80.11 






MC Simulation* [gj 


2.04 


2.07 


0.014 


0.15 


80.40 




bee 


Present result 


2.04 


2.07 


0.014 


0.17 


78.98 






MC Simulation* [2J,|2J 


2.03 


2.05 


0.010 


0.18 


78.40 


6 


fee 


Present result 


2.32 


2.35 


0.012 


0.12 


103.7 






MWDA-static reference [12| 


2.33 


2.35 


0.007 


0.17 








MC Simulation* [25j 
MC Simulation [2jj 


2.34 


2.37 


0.012 


0.15 


104.5 






2.32 


2.35 


0.012 




103.0 




bee 


Present result 


2.30 


2.33 


0.011 


0.16 


101.22 






MC Simulation* [2J 


2.32 


2.35 


0.011 


0.17 


103.6 






MSMC technique [23 
MC Simulation [26| 


2.30 


2.32 


0.011 




100.1 






2.30 


2.33 


0.012 




100.0 






MC Simulation* [2J,|2gl 


2.29 


2.31 


0.009 


0.18 


99.34 


4 


fee 


Present result 


5.60 


5.63 


0.008 


0.12 


565.6 






MWDA-static reference [12j 


5.22 


5.26 


0.008 


0.13 








MC Simulation [2J 


5.68 


5.71 


0.005 


0.17 


637.0 




bee 


Present result 


5.57 


5.61 


0.007 


0.16 


561.2 






MWDA-static reference [12| 


5.05 


5.09 


0.008 


0.18 








MC Simulation [2j| 


5.73 


5.75 


0.004 


0.18 


648.0 



NOTE-* indicates values obtained from interpolation of the tabulated values. 
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FIG. 1: Plots of &°'(r), ° a „ and — %^r-^ vs r for n = 6 and 7 = 2.30 which is close to the 



dp 



rV 



freezing point. The distance r is in unit of ag = I j 3 — I 
quantities for r > 1. 



1/3 



Insets show magnified values of respective 
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FIG. 2: Plot of t(r) vs r for n = 6 at 7 = 2.32 and n = 4 at 7 = 5.60. The distance r is in unit of 

/ o \l/3 

a ° = ( Ino I ■ r ^ ne dashed curve represents values for n = 4, 7; = 5.60 and full curve for n = 6, 
7* = 2.32. 




FIG. 3: Plot of —g vs r for n = 6 at 7 = 2.32 and n = 4 at 7 = 5.60. Other notations are same 
as in Fig. 2. 
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FIG. 4: Plot of s(r) vs r for n = 6 at 7 = 2.32 and n = 4 at 7 = 5.60. Other notations are same 
as in Fig. 2. 
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FIG. 5: Comparison of values of c, ' (r) as a function of r for a G vector of first set of fee and bec 

/ , \i/3 
lattices for n = 6, 7 S = 2.32, af cc = 32 and at, C c = 18. The distance r is in unit of ao = ( j^- j 

and fi = e~ > 4a . The dashed curve represents values of fee structure while full curve of bec 

structure. 
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(G 1) 
FIG. 6: Comparison of values of c\ ' (r) as a function of r for a G vector of second set of fee 

(dashed curve) and bee (full curve) lattices for n = 6, ^ s = 2.32, qj cc = 32 and a\, cc = 18. Other 

notations are same as in Fig. 5. 
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FIG. 7: Comparison of values of c ' (r) as a function of r for a G vector of the first six sets of 

/ s \i/3 
fee and bec lattices. The distance r is in unit of ao = ( -Sr- \ 
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FIG. 8: Comparison of values of c\ ),',(r) as a function of r for a G vector of the first set of fcc 



and bcc lattices for n = 6 at 7 S 



"linl'm' 

2.32, (Xf c 



32 and ab cc = 18. The distance r is in unit of 



«o 



I 4vrp ) 



\V3 



There is two sets of values for bcc lattice; one for G vectors lying in x-y plane and 



the other for the rest of G vectors of the first set. There is only one set of values for fcc lattice. 
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FIG. 9: Comparison of values of c ; m ;w (?") as a function of r for a G vector of the second set of fcc 
and bcc lattices. Other notations are same as in Fig. 8, except that there is now two sets of values 
shown by dashed and dotted curves for fcc lattice (see text) 



37 



3a = ka = 4.25; e G = 54.74 , 4> G = 45 



Ga = ka = 4.91 ; e e = , ty G = 90 




FIG. 10: Comparison of values (shown using color codes given on right hand side of each figure) of 
c (G,i) ^ g k ^ ^^ anc j C (G,2) ^ q^^ fc -j ag a f unc tion of cosdk (plotted on x-axis) and coscfik (plotted on 
y-axis) for Gia Q = ka = 4.25, 6 Gl = 54.7°, <j) Gl = 45° (a and c) and G 2 a = ka = 4.91, 6 G2 = 0°, 
4>c 2 = 90° (b and d) for a fee lattice for n = 6, 7 S = 2.32, aj cc = 32. 
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Ga = ka = 4.37; e G = 90 , <t> G = 45 
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(b) Ga = ka = 6.19;9 G = 90, <fe = 
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FIG. 11: Comparison of values (given in a color code) of & ' '(k,9k,<frk) an d c' ' 2 >(k,9k,4>k) as 
a function of cos6k and coscfrk for G\ao = kao = 4.37, 9q\ = 90°, (f>G 1 = 45° (a and c) and 
G2CL0 = kao = 6.19, 6g 2 = 90°, 4>q 2 = 0° (b and d) for a bcc lattice for n = 6, 7 S = 2.32, a^cc = 18. 
Other notation are same as in Fig. [10] 
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FIG. 12: Comparison of equilibrium phase diagram of — vs "fi found from simulation results and 



n 



from our theory. In inset the fluid-fee and fluid-bec transition lines are plated at a magnified scale 

and the fluid-bcc-fcc triple point is found at — = 0.158. 
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FIG. 13: Values of Cg (q, q, q) as a function of cosO (geometry is shown schematically in the figure) 
for qiao = 5200 = 4.3 and for potentials n = 12, 7/ = 1.17 (dashed curve) and n = 6, 7/ = 2.30 (full 
curve) 
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FIG. 14: Values of c 3 (q,q,q) vs qao (equilateral triangles). The dashed curve represents the 
values for n = 12, 72 = 1.17 and full curve for n = 6, 72 = 2.30. Inset shows values for qao > 4.0 on 
magnified scale. 
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FIG. 15: Values of c\ (<?i, (?2> I3) (shown using a color code shown on right hand side of each 
figure) as a function of <j) q2 and <f) q3 for q 1 = q 2 = q ma x with q ma xao 
(a) 9 q2 = 45°, ^ 3 = 45°, (b) 6 q2 = 45°, 8 qs = 90°, (c) ^ 2 = 90°, 9 q3 -45°, and (d) 9 q2 = 90°, 
(9 93 = 90° 



43 




360 
315 
270 

225 



• •- 



□ 



(b) 
















i 


i 













: 














































(d) 






i 


1 ' 




i 


l 


1 








-^ 




-\--_z 




































































i 
























i 









05 
0.4 



0.1 




0.5 





U: 



135 180 225 270 315 360 



90 135 180 225 270 315 360 



FIG. 16: Same as in Fig. [T5l except 6 qi = 90° and 



J <n 
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